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We  use  a  recently  developed  ab  initio  approach  to  calculate  the  lattice  thermal  conductivities  of  compound 
semiconductors.  An  exact  numerical  solution  of  the  phonon  Boltzmann  transport  equation  is  implemented,  which 
uses  harmonic  and  anharmonic  interatomic  force  constants  determined  from  density  functional  theory  as  inputs. 
We  discuss  the  method  for  calculating  the  anharmonic  interatomic  force  constants  in  some  detail,  and  we  describe 
their  role  in  providing  accurate  thermal  conductivities  in  a  range  of  systems.  This  first-principles  approach  obtains 
good  agreement  with  experimental  results  for  well-characterized  systems  (Si,  Ge,  and  GaAs).  We  determine  the 
intrinsic  upper  bound  to  the  thermal  conductivities  of  cubic  aluminum- V,  gallium- V,  and  indium-V  compounds  as 
limited  by  anharmonic  phonon  scattering.  The  effects  of  phonon-isotope  scattering  on  the  thermal  conductivities 
are  examined  in  these  materials  and  compared  to  available  experimental  data.  We  also  obtain  the  lattice  thermal 
conductivities  of  other  technologically  important  materials,  AIN  and  SiC.  For  most  materials,  good  agreement 
with  the  experimental  lattice  thermal  conductivities  for  naturally  occurring  isotopic  compositions  is  found.  We 
show  that  the  overall  frequency  scale  of  the  acoustic  phonons  and  the  size  of  the  gap  between  acoustic  and  optic 
phonons  play  important  roles  in  determining  the  lattice  thermal  conductivity  of  each  system.  The  first-principles 
approach  used  here  can  provide  quantitative  predictions  of  thermal  transport  in  a  wide  range  of  systems. 
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I.  INTRODUCTION 

Semiconductor  materials  have  played  an  important  role 
in  the  advancement  of  new  technologies  and  devices.  With 
the  growing  demand  for  more  efficient,  lower-cost,  reduced- 
size  electronic  devices,  consideration  and  manipulation  of 
the  constituent  material  properties  is  essential.  As  devices 
become  smaller  and  faster,  thermal  management  has  become 
an  increasingly  important  issue.  High  thermal  conductivities 
are  needed  for  passive  heat  spreading,  and  low  thermal  conduc¬ 
tivities  play  an  important  role  in  high  efficiency  thermoelectric 
materials.  In  semiconductors,  around  room  temperature  and 
higher,  lattice  vibrations  (phonons)  carry  the  majority  of  heat 
while  electronic  contributions  to  the  thermal  conductivity 
often  are  negligible.  The  lattice  thermal  conductivity  kl  is 
limited  mainly  by  intrinsic  phonon-phonon  scattering  which 
arises  from  crystal  anharmonicity  and  from  extrinsic  scattering 
processes  such  as  point  defects  and  boundaries.1 

To  understand  rigorously  how  material  properties  deter¬ 
mine  tty  and  to  devise  mechanisms  by  which  to  manipulate 
thermal  transport,  accurate  representation  of  the  intrinsic 
phonon-phonon  scattering  is  important.  Until  recent  years, 
full  microscopic  descriptions  of  this  scattering  have  been 
unavailable,  and  many  theories  of  kl  resorted  to  simple 
models  involving  a  number  of  ad  hoc  approximations.  Among 
them,  the  Debye  approximation  for  phonon  dispersions  was 
often  used,  neglecting  dispersion  in  the  acoustic  branches 
and  ignoring  optic  phonons  altogether,  and  mode  averaged 
Griineisen  constants  were  often  used  to  estimate  the  intrinsic 
phonon  scattering  rates. 1-5  Such  approximations  are  of  ques¬ 
tionable  validity,  and  because  such  models  are  typically  fit  to 
experimental  data,  they  lack  predictive  power.  Finally,  in  cases 
where  the  extrinsic  scattering  is  not  well  known,  these  models 
can  misrepresent  the  intrinsic  phonon-phonon  scattering  when 
obtaining  parameters  from  experiment. 


Recently,  first-principles  approaches  based  on  the  solution 
of  the  phonon  Boltzmann  transport  equation  (BTE)6  have  been 
developed,  which  provide  accurate  results  for  the  intrinsic 
phonon-phonon  scattering  and  kl  for  a  number  of  bulk 
semiconductor  systems,  and  have  demonstrated  good  agree¬ 
ment  with  measured  data  using  no  adjustable  parameters. 7-1 7 
In  these  approaches  interatomic  force  constants  (IFCs)  that 
are  required  to  calculate  phonon  frequencies  and  scattering 
rates  are  generated  using  density  functional  theory1819  (DFT) 
and/or  density  functional  perturbation  theory  (DFPT).20  This 
approach  has  the  advantage  of  being  predictive  and  can  readily 
and  accurately  be  applied  to  many  materials  in  contrast  to 
empirical  interatomic  potentials.  We  note  that  other  methods 
such  as  molecular  dynamics  within  the  Green-Kubo  formalism 
have  also  been  developed  to  calculate  kl  using  IFCs  from 
DFT.21  In  some  sense,  the  BTE  and  molecular  dynamics 
approaches  are  complementary:  The  BTE  directly  calculates 
quantum  mechanical  scattering  rates,  but  it  includes  only 
lowest  order  phonon-phonon  interactions  and  so  becomes 
less  accurate  at  very  high  temperatures  where  higher  order 
processes  can  become  important.  The  molecular  dynamics 
approach  includes  anharmonicity  to  all  orders,  but  it  is  based 
on  a  classical  picture  and  so  less  valid  at  lower  temperatures. 
Accurate  results  for  kl  with  molecular  dynamics  methods 
have  been  obtained  and  depend  on  system  size  and  sampling 
methods,  which  can  become  computationally  costly.  Until  now 
first-principles  methods  have  been  applied  to  relatively  few 
systems,  including  silicon  (Si),  germanium  (Ge),  and  their 
alloys,7,8  diamond, 9-11  GaAs,12  magnesium  oxide  (MgO),13 
half-Heusler  compounds, 14  lead  selenide  (PbSe),  lead  telluride 
(PbTe),  and  their  alloys, 15  gallium  nitride  (GaN),16  magnesium 
silicide  (Mg2Si),  magnesium  stannide  (MgjSn),  and  their 
alloys.17  In  this  paper  we  apply  a  first-principles  BTE  approach 
to  determine  kl  in  a  range  of  technologically  important 
compound  semiconductors. 
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In  Sec.  II  we  briefly  describe  the  BTE  thermal  transport 
theory  and  discuss  details  of  our  real-space  DFT  calculations 
for  determining  the  anharmonic  IFCs  that  govern  phonon- 
phonon  scattering.  The  calculated  k\  for  the  test  cases,  Si, 
Ge,  and  GaAs,  are  presented  in  Sec.  Ill  with  discussion  of  the 
effects  of  supercell  size,  enforcement  of  symmetry  conditions, 
and  interaction  radii  on  the  anharmonic  IFCs  and  at l.  Results 
and  discussion  are  given  for  aluminum- V,  gallium- V,  and 
indium-V  cubic  compounds  in  Sec.  IV  and  for  cubic  silicon 
carbide  (3C-SiC)  and  wurtzite  aluminum  nitride  (AIN)  in 
Sec.  V.  A  summary  and  conclusions  are  given  in  Sec.  VI.  We 
compare  the  calculated  phonon  dispersions  for  the  materials 
studied  in  this  paper  with  available  experimental  data  in  the 
Appendix. 

II.  THERMAL  TRANSPORT  AND  ANHARMONIC  IFCs 

Quantitative  understanding  of  atl  requires  an  accurate 
microscopic  description  of  the  intrinsic  anharmonic  phonon- 
phonon  scattering.  Here  we  briefly  discuss  the  phonon 
Boltzmann  transport  equation  approach  used  to  determine 
rigorously  k\  .  Details  of  this  approach  can  be  found  in 
Refs.  7,9-11,16,  and  22-24.  The  thermal  conductivity  tensor 
is 

Kgfi  =  -^7  (dn°x/dT)na>x.VxaVx.pr>lCl,  (1) 

where  V  is  the  crystal  volume,  a  and  /3  are  Cartesian  directions, 
is  the  Bose  distribution  function,  T;„  is  the  phonon 
lifetime,  cox  is  the  phonon  frequency,  and  Vxa  =  da>x/dqa  is 
the  component  of  the  phonon  velocity  in  the  ath  direction. 
Here,  X  =  ( q,j )  designates  a  phonon  with  wave  vector  q  in 
branch  j .  For  the  systems  examined  here  Kap  is  diagonal  and 
can  be  described  by  a  scalar  for  cubic  structures  and  by  in- 
and  out-of-plane  components  for  wurtzite  structures.25  The 
phonon  frequencies  are  determined  by  diagonalization  of  the 
dynamical  matrix  (see  the  Appendix). 

We  consider  k\  to  be  limited  by  intrinsic  three-phonon 
scattering  and  by  extrinsic  point  defect  scattering  by  isotopic 
impurities.  For  high  quality  semiconductors  these  are  the 
dominant  scattering  mechanisms  for  phonons  above  a  few 
tens  of  kelvin.1  The  intrinsic  anharmonic  scattering  rates 
l/tfh  (Ref.  16)  are  determined  from  scattering  processes 
involving  three  phonons  that  satisfy  conservation  of  energy 
and  momentum:  cox  ±  my  =  a>x"  and  q  ±  q'  =  q"  +  K .  The 
±  signs  correspond  to  the  two  types  of  possible  three-phonon 
processes,  and  K  is  a  reciprocal  lattice  vector,  which  is 
zero  for  normal  (N)  processes  and  nonzero  for  umklapp  (U) 
processes.1  The  isotope  scattering  rates  I  / r-ls0  are  determined 
from  perturbation  theory.26,27  The  phonon  scattering  time 
within  the  single  mode  relaxation  time  approximation  (RTA) 

rRTA  _  [i/ranh_|_  i/T|so]-t  can  determined  fr0m  these 

scattering  rates  alone,  however,  the  full  solution  of  the  BTE 
is  required  to  determine  the  actual  phonon  transport  time 
r^. 7,9-11,16,22-24,28,29  Using  r^TA  in  the  calculation  of  the 
thermal  conductivity  krta.  the  N  phonon-phonon  scattering 
processes  are  incorrectly  treated  as  resistive.1  The  full  BTE 
solution  corrects  for  this  and  therefore  gives  a  larger  K\  .  The 
difference  between  K\  from  the  full  solution  and  arta  gives 
a  measure  of  the  importance  of  the  N  scattering  relative  to 


U  scattering.  In  most  of  the  materials  considered  here,  the 
U  scattering  is  relatively  strong  around  room  temperature, 
so  krta  does  not  differ  much  from  the  full  solution  to  the 
BTE.  Similar  behavior  has  been  previously  noted  in  Si  and 
Ge.7,30  For  systems  with  very  strong  N  scattering  relative 
to  U  scattering,  such  as  in  diamond, 9-11  graphene,31'2  and 
carbon  nanotubes,22,33  the  full  solution  to  the  BTE  is  required 
to  accurately  determine  A|  . 

The  BTE  formalism  requires  harmonic  IFCs,  <bap(lK,l'  tc'), 
as  inputs  to  determine  the  phonon  frequencies,  eigenvectors, 
and  velocities.  The  anharmonic  IFCs,  <&apy(lK,l' k' ,1" k"), 
are  required  to  determine  the  phonon-phonon  scattering 
rates. 7,9-1 1,16,22-24  All  IFCs  were  calculated  using  norm- 
conserving  pseudopotentials  in  the  local  density  approx¬ 
imation  (LDA)  with  the  plane-wave  QUANTUM  ESPRESSO 
package.34,35  The  ground-state  configuration  of  each  system 
was  determined  by  adjusting  the  structural  lattice  constants 
(a  for  cubic,  a  and  c  for  wurtzite)  to  find  the  minimum  energy.36 
The  harmonic  IFCs  were  calculated  using  standard  DFPT. 
Typically,  convergence  was  achieved  with  80  Ryd  plane- 
wave  cutoffs  and  6x6x6  A:-point  meshes  for  electronic  and 
phonon  calculations. 

Calculation  of  the  anharmonic  IFCs  is  not  part  of  typical 
electronic  structure  packages.  The  anharmonic  IFCs  used  in 
this  work  were  calculated  using  a  real-space  “frozen  phonon” 
approach  within  DFT.16,37  We  define  a  nearest  neighbor  cutoff 
radius  for  which  interactions  between  atoms  separated  by  a 
distance  larger  than  the  cutoff  are  taken  to  be  zero.  A  large  num¬ 
ber  of  anharmonic  three-phonon  elements,  (0 k,V k'  ,1" k"), 

need  to  be  determined,  which  can  be  computationally  costly. 
Fortunately,  symmetry  considerations  establish  relationships 
between  IFCs  and  limit  the  number  of  independent  elements 
to  be  calculated.37-39  We  use  symmetry  conditions  to  identify 
a  set  of  independent  IFCs,  0, ,  from  which  all  anharmonic  IFCs 
can  be  linearly  constructed: 

k' ,1" k")  =  Y  4>iAi(lic,r k' ,1" k").  (2) 

i 

The  sum  in  Eq.  (2)  runs  over  all  independent  anharmonic 
IFCs  and  Ai(hc,l'Kr ,1"k")  are  constants  determined  by  the 
symmetry  operations  that  establish  the  relationship  of  each 
( pi  with  each  anharmonic  IFC.  For  cubic  crystals  there  is  a 
simple  one-to-one  correspondence  between  anharmonic  IFCs 
giving  Ai(lK,l' k' d" k")  —  0  for  all  but  one  element  for  which 
A, (/a:, /V, /"&■")  =  ±1.  For  structures  with  less  symmetry, 
such  as  hexagonal  lattices,  there  is  not  necessarily  a  one-to-one 
correspondence  between  IFCs  and  the  Ai(hc,l'ic'  are  not 
simple  integers. 

For  each  independent  IFC  element,  pairs  of  atoms  (one 
from  the  unit  cell  and  one  of  its  neighbors  depending  on 
the  element)  were  systematically  perturbed  from  equilibrium 
by  a  small  distance  within  a  large  ground-state  supercell. 
These  supercells  must  be  sufficiently  large  so  that  interactions 
with  perturbed  periodic  supercell  images  do  not  contribute 
significantly  to  the  calculations.  The  resulting  Hellmann- 
Feynman  forces34  on  all  the  atoms  were  calculated  via  Y -point 
self-consistent  calculations  for  four  different  perturbations 
to  obtain  numerical  derivatives  of  the  forces  and  thus  the 
anharmonic  IFC  element.  Since  the  forces  for  all  of  the 
atoms  are  calculated  for  a  given  perturbation,  this  method 
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overdetermines  the  independent  IFC  elements.  We  find  that 
numerical  differences  between  calculated  elements  that  should 
be  the  same  by  symmetry  are  typically  small  and  subsequently 
use  their  average. 

Before  applying  symmetry  operations  to  the  independent 
IFC  elements  to  determine  the  dependent  IFCs,  other  important 
symmetry  properties  need  to  be  considered.  In  particular,  the 
truncation  of  anharmonic  interactions  to  a  finite  number  of 
neighbors  breaks  translational  invariance  (TI)  of  the  crystal. 
Translational  invariance  can  play  an  important  role  in  deter¬ 
mining  the  scattering  properties  of  long-wavelength  acoustic 
phonons  and  thus  in  determining  k\  .  The  TI  conditions  can  be 
written  as17'38 

tjAij  =  Bt  =  0.  (3) 

V'k"  j 

Due  to  derivative  permutation  symmetry  there  is  also  a  set  of 
equations  for  the  sum  over  I'k' .  Each  i  represents  a  separate  TI 
equation  for  each  combination  of  variables  that  is  not  summed 
over.  B,  is  the  sum  for  each  TI  condition  that  should  be  zero  for 
perfect  translational  invariance,  and  we  use  the  more  compact 
notation  A,j  =  J2t'K'  or  i»K«  Ay-((k,/V,/'V0  to  write  Eq.  (3)  in 
terms  of  the  independent  IFCs,  </);.  We  employed  and  tested 
three  different  procedures  for  enforcing  the  TI  conditions:  (i) 
simple  acoustic  sum  rule  (ASR),  (ii)  y2  minimization,  and 
(iii)  Lagrange  multipliers.  For  the  ASR  method,  a  subset  of 
anharmonic  IFCs  is  defined  in  terms  of  the  others, 

«W0*,/V,0k)  =  -  J2  (4) 

1"k"^0k 

as  well  as  permutations.  Using  the  ASR  method  has  the  advan¬ 
tage  of  being  simple  and  guarantees  that  the  TI  conditions  are 
met.  However,  it  does  not  enforce  derivative  permutation  and 
point  group  symmetries  of  the  system.  For  the  second  method 
we  use  a  numerical  algorithm4"  to  minimize  the  function 
X2  —  B2  by  adding  a  small  compensation  to  each  </), .  This 
procedure  and  the  Lagrange  multipliers  method  (described  in 
Refs.  17  and  41)  change  only  the  independent  IFCs  and  thus 
maintain  all  symmetry  properties  of  the  crystal  while  enforcing 
the  TI  conditions.  Given  the  set  of  independent  IFCs,  the 
symmetry  operations  can  then  be  used  to  determine  all  other 
anharmonic  IFCs. 

m.  TEST  CASES  (Si,  Ge,  AND  GaAs) 

Here  we  present  calculated  results  for  atl  of  Si,  Ge, 
and  GaAs  and  quantify  our  discussion  of  various  aspects 
of  the  anharmonic  calculations  in  terms  of  their  effects 
on  K\  .  These  materials  have  been  extensively  studied  and 
reliable  measurements  of  K[  for  both  isotopically  enriched 
and  naturally  occurring  isotope  concentrations  are  available. 

Figure  1  shows  the  calculated  atl  versus  temperature  for  Si. 
Solid  curves  are  calculated  kl  for  isotopically  pure  Si,  /cpure, 
and  dashed  curves  are  calculated  atl  with  naturally  occurring 
Si  isotope  concentrations,  Arnaturai.  Circles42  are  experimental 
atl  for  isotopically  enriched  Si  and  triangles42  and  squares43,44 
are  experimental  ^natural-  We  note  that  the  atl  data  from  Ref.  42 
agree  well  with  the  data  from  Ref.  45  (not  shown)  taken 
from  three  different  laboratories.  The  black  and  red  curves 
are  for  a:l  determined  with  anharmonic  IFCs  calculated  in 


PHYSICAL  REVIEW  B  87,  165201  (2013) 


temperature  (K) 


FIG.  1.  (Color  online)  Calculated  atl  vs  temperature  for  Si  with 
experimental  data  for  isotopically  enriched  Si  [circles  (Ref.  42)]  and 
for  naturally  occurring  Si  concentrations  [triangles  (Ref.  42)  and 
squares  (Refs.  43  and  44)].  Solid  curves  correspond  to  calculated 
/Cpure  and  dashed  curves  correspond  to  Arnatliral.  was  determined  with 
anharmonic  IFCs  calculated  in  a  64  atom  supercell  (black  curves)  and 
a  216  atom  supercell  (red  curves).  The  calculated  Arnaturai  for  a  64  atom 
supercell  without  enforcing  TI  is  given  by  the  green  dashed  curve. 

a  64  atom  supercell  and  a  216  atom  supercell,  respectively. 
All  calculations  include  anharmonic  IFCs  out  to  third  nearest 
neighbors  and  the  TI  conditions  enforced  using  the  y2 
minimization  procedure.  The  green  dashed  curve  (lowest)  is 
^natural  determined  with  anharmonic  IFCs  calculated  in  a  64 
atom  supercell  but  without  enforcing  the  symmetry  conditions. 
The  red  curves  (216  atom  supercell)  show  excellent  agreement 
with  the  /cl  data  from  Ref.  42  over  a  wide  temperature  range 
without  use  of  adjustable  parameters.  We  also  calculated 
kl  for  isotopically' purified  Si  (99.983%  28Si,  0.014%  29Si, 
and  0.003%  30Si)  corresponding  to  experiment42  and  find 
negligible  difference  with  Arpure  over  the  entire  temperature 
range.  The  black  curves  (64  atom  supercell)  demonstrate  good 
agreement  with  the  red  curves  even  though  a  much  smaller 
supercell  was  used  to  calculate  the  anharmonic  IFCs.  The  green 
curve  for  A-naturai  without  enforcing  the  TI  conditions  on  the 
anharmonic  IFCs  is  significantly  lower  than  the  red  and  black 
curves  and  experiment.  The  absence  of  TI  causes  an  artificial 
enhancement  of  low  frequency  scattering  of  heat-carrying 
acoustic  phonons  and  thus  reduces  atl- 

Figure  2  shows  the  calculated  Arpure  (solid  curves)  and 
^natural  (dashed  curves)  for  Ge  with  experimental  values  for 
isotopically  purified  Ge  (circles)  and  with  naturally  occur¬ 
ring  isotopic  abundances  (triangles)  46  The  red,  black,  and 
purple  curves  correspond  to  atl  determined  with  translational 
invariance  enforced  via  the  /2  minimization  procedure,  the 
Lagrange  multiplier  method,  and  the  simple  ASR  method, 
respectively.  The  green  dashed  curve  corresponds  to  Arnaturai 
without  symmetry  conditions  enforced  on  the  anharmonic 
IFCs.  As  with  Si,  the  Arnaturai  of  Ge  as  calculated  without 
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FIG.  2.  (Color  online)  Calculated  a:l  vs  temperature  for  Ge  with 
experimental  data  for  isotopically  enriched  Ge  (circles)  and  for 
naturally  occurring  Ge  concentrations  (triangles)  (Ref.  46).  kl  was 
determined  with  TI  enforced  via  the  x2  minimization  procedure 
(red  curves),  the  Lagrange  multiplier  method  (black  curves),  and  the 
simple  ASR  method  (purple  curves).  fcpure  is  given  by  solid  curves  and 
^natural  is  given  by  dashed  curves.  The  green  dashed  curve  corresponds 
to  /r natural  without  TI  conditions  enforced. 

enforcing  the  TI  conditions  is  much  less  than  experiment  and 
the  other  calculations.  The  calculated  /rnatuiai  using  the  simple 
ASR  (purple  dashed  curve)  to  enforce  TI,  which  results  in 
violation  of  other  crystal  symmetries,  also  falls  significantly 
below  both  experiment  and  the  corresponding  black  and  red 
curves.  The  k\  calculated  using  the  x2  minimization  procedure 
(red  curves)  and  the  Lagrange  multipliers  method  (black 
curves)  are  in  good  agreement  with  each  other  and  with 
experiment.  Figures  1  and  2  demonstrate  that  TI  plays  an 
important  role  in  determining  k l  and  show  that  as  long  as  TI 
and  the  point  group  and  permutation  symmetries  are  satisfied, 
the  method  for  TI  enforcement  is  not  a  critical  issue. 

We  test  the  effects  of  the  nearest  neighbor  cutoff  radius 
on  for  GaAs  in  Fig.  3.  Since  GaAs  is  a  polar  compound 
with  long  range  Coulomb  interactions,  the  calculation  of  k\ 
in  GaAs  may  be  more  sensitive  to  the  cutoff  radius  than  Si 
or  Ge.  Figure  3  shows  the  calculated  /cnaturai  for  GaAs  with 
different  nearest  neighbor  cutoff  radii  along  with  experimental 
data  for  /cnaturai  (green  triangles47  and  black  triangles48)  and 
kl  for  isotopically  purified  GaAs  (green  circles47  and  black 
circles48).  The  calculated  curves  included  second  (green),  third 
(solid  black),  fourth  (red),  and  fifth  (blue)  nearest  neighbors 
for  the  anharmonic  IFCs.  The  second,  third,  fourth,  and  fifth 
nearest  neighbor  cutoff  radii  were  0.75a,  0.85 a,  1.05a,  and 
1.1a,  respectively,  which  included  interactions  of  unit  cell 
atoms  with  16,  28,  34,  and  46  surrounding  atoms,  respectively. 
Here  a  is  the  lattice  constant  (see  Table  I).  We  found  that 
including  only  first  nearest  neighbors  required  very  large 
changes  to  the  anharmonic  IFCs  in  order  to  enforce  the  TI 
conditions  and  thus  was  not  included  here.  We  have  also  used 


FIG.  3.  (Color  online)  Calculated  /fnaturai  vs  temperature  for  GaAs 
with  experimental  data  for  isotopically  enriched  GaAs  [green  circles 
(Ref.  47)  and  black  circles  (Ref.  48)]  and  for  naturally  occurring 
Ga  concentrations  [green  triangles  (Ref.  47)  and  black  triangles 
(Ref.  48)].  The  calculated  curves  included  second  (green),  third  (solid 
black),  fourth  (red),  and  fifth  (blue)  nearest  neighbors  for  determining 
the  anharmonic  IFCs.  The  dashed  black  curve  gives  /cnatural  for  GaAs 
using  IFCs  determined  from  DFPT  to  seventh  nearest  neighbors  and 
with  long  range  Coulomb  interactions  included  (described  in  text). 

a  reciprocal-space  DFPT  approach  (black  dashed  curve)  to 
calculate  the  anharmonic  IFCs  and  to  calculate  ^natural,  which 
extends  the  interactions  to  seventh  nearest  neighbors  and 
includes  long  range  Coulomb  interactions.7,9-11,49,50  As  can 
be  seen  in  Fig.  3,  the  calculated  ^natural  is  fairly  insensitive  to 
the  nearest  neighbor  cutoff  radius.  The  inset  to  Fig.  3  gives  an 
expanded  view  of  the  curves  around  room  temperature.  With 
an  increasing  number  of  nearest  neighbors,  the  magnitudes  of 
the  anharmonic  IFCs  obtained  from  both  real-space  and  DFPT 
approaches  involving  distant  atoms  become  smaller,  which 
provides  more  potential  for  accuracy  errors.  Nevertheless,  the 
variation  in  the  room  temperature  GaAs  kl  values  shown  for 
the  different  cases  in  the  inset  to  Fig.  3  is  less  than  4%. 
We  find  that  the  DFPT  approach  gives  a  similar  variation  of 
atl  for  different  numbers  of  nearest  neighbors  (not  shown). 
This  suggests  a  robustness  of  the  first-principles  real-space 
approach  for  k\  presented  here.  All  of  the  theoretical  curves  lie 
above  the  experimental  data  from  Ref.  48  and  agree  somewhat 
better  with  the  data  from  Ref.  47.  We  also  note  that  K\  of  GaAs 
does  not  have  much  enhancement  with  isotopic  purification, 
~5%  at  room  temperature. 

We  have  shown  through  these  test  cases  (Si,  Ge,  and  GaAs) 
that  calculations  of  k\  are  fairly  insensitive  to  reasonable 
choices  of  supercell  size  and  nearest  neighbor  cutoff  distance. 
We  have  also  demonstrated  that  it  is  important  to  enforce 
translational  invariance  conditions  while  maintaining  the  point 
group  and  derivative  permutation  symmetries  to  accurately 
determine  K\  . 
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IV.  III-V  CUBIC  SEMICONDUCTORS 

In  this  section  we  apply  the  first-principles  approach  for 
atl  discussed  in  Sec.  II  and  tested  on  well-studied  Si,  Ge,  and 
GaAs  to  the  III-V  cubic  compound  materials.  For  all  cases,  we 
used  216  atom  supercells  with  an  interaction  cutoff  to  include 
third  nearest  neighbors  for  calculating  anharmonic  IFCs.  We 
enforced  translational  invariance  via  the  x2  minimization  pro¬ 
cedure.  The  lattice  constants  used  to  calculate  the  harmonic  and 
anharmonic  IFCs  were  determined  by  energy  minimization 
except  for  InP,  InAs,  and  InSb  (discussed  below).  For  each 
material  in  this  paper,  the  lattice  constants  and  mass  variance 
parameters26  used  in  the  calculations  of  al  and  the  calculated 
KVme,  tr natural.  and  percent  isotope  effect  (percent  increase  in 
A:pUre  compared  to  /(natural )  p  —  100  x  (/Cpure/zc natural  -  1)  are  in 
Table  I.  Also  in  Table  I  is  the  percent  increase  to  the  room 
temperature  Apure  from  the  full  solution  to  the  BTE  compared 
to  the  RTA  solution  S  —  100  x  (/(pure/ arta  —  1)-  For  small  S, 
U  phonon-phonon  scattering  is  large  and  the  RTA  solution 
to  the  BTE  is  a  good  approximation  for  al-  For  large  .S'.  N 
phonon-phonon  scattering  processes  are  important  and  the 
fully  iterated  solution  is  required  to  accurately  give  al- 

The  calculated  Apure  for  the  different  materials  is  limited 
by  only  three-phonon  scattering  and  represents  the  intrinsic 


upper  bound  to  al-  To  compare  better  with  experiment  we  also 
include  isotopic  impurity  scattering  to  determine  AnatUi-ai  for 
systems  composed  of  elements  with  differing  natural  isotopic 
abundances.  We  note  that  experimental  samples  may  contain 
defects  (dislocations,  substitutions,  other  impurities,  etc.). 
Further,  at  low  temperatures,  sample  size  and  grain  boundaries 
can  play  a  significant  role  in  determining  al-  These  scattering 
mechanisms  are  not  considered  in  this  work,  but  they  can 
significantly  lower  al  below  the  values  obtained  here. 

A.  Aluminum-V  compounds 

Here  we  examine  al  of  A1P,  AlAs,  and  AlSb.  Figure  4 
shows  calculated  al  versus  temperature  for  isotopically  pure 
A1P  (solid  green  curve),  AlAs  (solid  red  curve),  and  AlSb 
(solid  black  curve).  Phonon-isotope  scattering  is  particularly 
important  in  AlSb  in  part  due  to  the  large  isotope  mixture  of  Sb 
(57.21%  121  Sb  and  42.79%  123 Sb)  while  the  isotope  mixtures 
for  Al,  P,  and  As  are  negligible.  Thus,  the  calculated  Anaturai 
is  shown  only  for  AlSb  by  the  dashed  black  curve  in  Fig.  4. 
The  circles  correspond  to  experimental  data  for  ArnatUi-ai  of  A1P 
(solid  green51)  and  AlSb  (solid  black51  and  open  black52).  An 
experimental  data  point  for  al  of  AlAs  at  7'  =  300  K  (Ref.  53) 
(not  shown)  lies  directly  on  the  solid  green  circle.  Optic 


TABLE  I.  Experimental  and  calculated  lattice  constants  determined  by  energy  minimization,  mass  variance  parameters,  Apllre,  Anatural,  and 
percent  isotope  effect  (P)  at  T  =  300  K  for  materials  considered  in  this  work.  Also  included  is  a  comparison  of  the  room  temperature 
Apure  determined  by  the  BTE  and  the  thermal  conductivity  determined  by  the  RTA,  arta,  for  each  material.  The  difference  is  given  by 
S  =  100  X  (Apure/ARTA  !)• 


^calc  (^exp)  (^) 

^cation  (<?anion)  (  X  10  ) 

Apure  (Win-1  K"1) 

Anatural  (Wm"1  KG1) 

P 

S 

Diamond" 

3.53  (3.57) 

0.75  (0.75) 

3450 

2290 

51 

50 

Si 

5.37  (5.43) 

2.01  (2.01) 

155 

144 

8 

2.4 

Ge 

5.61  (5.65) 

5.87  (5.87) 

74 

60 

23 

6.9 

3C-SiC 

4.34  (4.36) 

2.01  (0.75) 

572 

479 

20 

6.7 

A1P 

5.40  (5.45) 

-  (-) 

90 

- 

- 

3.0 

AlAs 

5.61  (5.66) 

-  (-) 

105 

- 

- 

0.5 

AlSb 

6.10(6.14) 

-  (0.66) 

118 

86 

36 

39 

c-GaNb 

4.42  (4.50) 

1.97  (-) 

362 

215 

68 

13 

GaP“ 

5.34  (5.45) 

1.97  (-) 

153 

131 

16 

4.8 

GaAsc 

5.55  (5.65) 

1.97  (-) 

56 

54 

4 

5.6 

GaSb“ 

6.00  (6.10) 

1.97  (0.66) 

48 

45 

6 

3.1 

InP1’ 

5.79  (5.87) 

0-12  ( — ) 

91 

89 

2 

2.4 

InAsb 

5.97  (6.06) 

0-12  ( — ) 

36 

36 

.5 

2.8 

InSbb 

6.39  (6.48) 

0.12(0.66) 

20 

20 

2 

3.9 

U)-GaNb-d 

3.13  (3.19) 
5.10(5.19)“ 

0.377  (0.377)f 

1.97  (-) 

401  (385/ 

242  (239)s 

66 

7.2 

uj-AIN 

3.05  (3.11) 

4.81  (4.98)“ 

0.387  (0.382/ 

-  (-) 

322  (303/ 

14 

“Calculated  in  Ref.  1 1 . 


bThese  calculated  lattice  constants  determined  by  energy  minimization  were  increased  to  better  fit  the  phonon  dispersions  when  calculating  al 
for  these  systems,  as  explained  in  the  text. 

“The  al  and  P  values  presented  here  slightly  differ  from  the  calculated  results  from  Ref.  16  which  used  ~0.2%  increase  to  the  lattice  constants 
to  account  for  zero-point  and  finite  temperature  atomic  motion. 
dCalculated  in  Ref.  16. 

“Calculated  and  experimental  lattice  constant  c  for  the  wurtzite  structure. 
fCalculated  and  experimental  internal  parameter  u  for  the  wurtzite  structure. 

EOut-of-plane  component  to  al  along  the  c  axis. 
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FIG.  4.  (Color  online)  Calculated  Arpure  vs  temperature  for  A1P 
(solid  green  curve),  AlAs  (solid  red  curve),  and  AlSb  (solid  black 
curve).  Experimental  data  for  kl  of  A1P  are  solid  green  circles 
(Ref.  51)  and  for  kl  of  AlSb  are  solid  (Ref.  51)  and  open  (Ref.  52) 
black  circles.  An  experimental  data  point  for  kl  of  AlAs  at  T  =  300  K 
(Ref.  53)  (not  shown)  lies  directly  on  the  solid  green  circle.  Calculated 
^natural  for  AlSb  is  given  by  the  dashed  black  curve.  The  acoustic-only 
contribution  to  /cnat urai  for  AlSb  is  given  by  the  dotted  black  curve.  The 
inset  gives  /cpure  (black  curve),  the  acoustic  contribution  to  >cpure  (red 
curve),  and  the  optic  contribution  to  Arpure  (green  curve)  vs  temperature 
for  AlSb. 

phonons  play  an  important  role  in  determining  k\  for  all  of 
the  materials  considered  in  this  work,  but  especially  so  when 
comparing  K\  of  the  aluminum-V  compounds.  We  note  that 
the  importance  of  optic  phonon  scattering  of  heat-carrying 
acoustic  phonons  has  been  previously  discussed  in  a  number 
of  papers,  including  Refs.  11,  15-17,  22,  28,  29,  51,  and  54. 
In  AlSb  the  large  mass  mismatch  of  A1  and  Sb  atoms  results 
in  a  substantial  frequency  gap  between  acoustic  and  optic 
branches  (see  Fig.  13  in  the  Appendix).  Important  phonon 
scattering  processes  of  the  type  acoustic  +  acoustic  -o-  optic 
( aao )  are  completely  forbidden  by  energy  conservation 
due  to  the  large  gap  in  AlSb.  Further,  energy  conservation 
forbids  ooo  scattering  in  all  the  materials  considered  here 
and  the  phase  space  for  aoo  scattering  in  AlSb  is  very  small 
due  to  the  frequency  gap.  These  restrictions  become  less 
severe  as  the  mass  differences  between  A1  and  anions  decrease 
and  the  frequency  gaps  decrease51'55  (see  Figs.  11-13  in  the 
Appendix). 

At  T  =100  K  the  optic  phonons  are  not  significantly 
thermally  populated,  and  so  they  do  not  provide  much 
scattering  resistance  to  the  lower  frequency  acoustic  phonons 
regardless  of  the  frequency  gap.  In  this  temperature  regime 
aaa  scattering  provides  the  dominant  thermal  resistance  in 
all  materials.  A1P  generally  has  the  weakest  aaa  scattering 
of  these  compounds  due  to  a  higher  frequency  scale  of  the 
acoustic  phonons,  which  decreases  phonon  populations  and 
limits  the  phase  space  for  this  scattering.  Further,  A1P  has 


the  highest  acoustic  velocities.  Thus,  at  100  K,  /rpure  for  A1P 
is  2.2  times  higher  than  /cpure  for  AlSb  and  is  the  highest  of 
the  three  materials.  With  increasing  temperature  anharmonic 
phonon  scattering  becomes  greater  and  a:l  decreases,  as  can 
be  seen  for  all  the  curves  in  Fig.  4.  The  /rpure  for  A1P  decreases 
faster  with  increasing  temperature  than  Arpure  for  AlAs  and 
AlSb  due  to  stronger  aao  scattering  in  A1P  as  optic  phonons 
are  increasingly  thermally  populated.  At  T  —  300  K,  A1P  has 
the  lowest  Arpure  of  the  three  materials  due  to  having  the  smallest 
frequency  gap  and  thus  stronger  acoustic-optic  scattering. 

AlSb  has  the  highest  k l  at  room  temperature  due  to  the 
severely  restricted  acoustic-optic  scattering  compared  to  the 
other  materials,  which  leads  to  both  increased  acoustic  phonon 
lifetimes  and  increased  optic  phonon  lifetimes.  Unlike  all  of  the 
other  materials  in  this  work,  the  phase  space  for  three-phonon 
scattering  of  optic  phonons  in  AlSb  is  so  restricted  that  the 
optic  modes  provide  significant  contributions  to  atl  despite 
having  lower  velocities  than  typical  acoustic  phonons.  Since 
aao  and  ooo  scattering  channels  are  completely  forbidden  by 
energy  conservation,  aoo  scattering  provides  the  only  intrinsic 
resistance  for  optic  phonons  in  AlSb.  The  phase  space  for  aoo 
scattering  is  sensitive  to  the  splitting  of  the  longitudinal  optic 
(LO)  branch  and  the  transverse  optic  (TO)  branch  which  arises 
from  long  range  Coulomb  interactions.  Increased  LO/TO 
splitting  leads  to  a  larger  phase  space  for  aoo  scattering  and 
reduces  optic  phonon  contributions  to  /cL,  while  having  little 
effect  on  the  acoustic  phonon  contributions. 

We  note  that  in  PbTe,15  Mg2Si,17  and  single-walled  carbon 
nanotubes,22’33  optic  phonon  contributions  to  K{  were  also 
found  to  be  significant.  However,  in  those  cases  the  large 
contributions  stemmed  from  optic  branches  that  were  relatively 
dispersive  and  resided  at  lower  frequencies  rather  than  from 
a  restricted  phase  space  for  scattering.  The  inset  to  Fig.  4 
shows  the  /rpure  versus  temperature  for  AlSb  (black  curve) 
with  the  contributions  to  /rpure  from  the  acoustic  phonons  (red 
curve)  and  the  optic  phonons  (green  curve).  For  T  >  250  K 
the  optic  modes  actually  provide  the  dominant  contribution 
to  K]  ,  As  the  temperature  decreases,  the  contributions  to  kl 
from  both  acoustic  and  optic  modes  increase  as  anharmonic 
scattering  gets  smaller,  but  the  acoustic  phonon  contribution  to 
atl  continues  to  increase  while  the  optic  phonon  contribution 
peaks  at  T  —  1 75  K  and  then  drops.  This  reflects  the  significant 
decrease  in  the  optic  phonon  population  with  decreasing 
temperature  due  to  their  higher  frequencies,  which  reduces 
their  contributions  to  aj  .  The  large  gap  in  AlSb  substantially 
weakens  the  U  scattering  compared  to  N  scattering  for  the  optic 
modes.  Thus,  Arpure  from  the  full  solution  to  the  BTE  is  40% 
higher  than  a:rta  (see  Table  I).  The  acoustic  contributions  to 
kl  of  AlSb  only  give  S  =  2%,  while  the  optic  modes  give  S  — 
1 16%.  In  contrast,  the  total  increase  is  only  S  =  3%  for  A1P. 

We  find  good  agreement  with  the  experimental  atl  data 
for  A1P  and  AlAs,  however,  the  calculated  Arnaturai  for  AlSb 
is  ~50%’  above  experiment  at  T  —  300  K.  The  experimental 
atl  data  was  taken  from  an  AlSb  sample  believed  to  have 
had  appreciable  amounts  of  impurities.52  We  note  that 
point  defects  can  be  particularly  effective  at  scattering 
high  frequency  optical  phonons.  The  dotted  black  curve  in 
Fig.  4  gives  the  acoustic-only  contribution  (assuming  that 
additional  defects  eliminate  the  optic  contributions  as  in 
typical  materials)  to  Arnaturai  of  AlSb,  which  is  in  excellent 
agreement  with  the  experimental  data. 
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FIG.  5.  (Color  online)  Calculated  /cL  vs  temperature  for  GaSb 
(black  curve),  GaAs  (green  curve),  GaP  (red  curves),  and  cubic  GaN 
(blue  curves).  The  solid  curves  give  /cpure  and  the  dashed  curves 
give  /cnatUrai  •  Experimental  data  is  shown  for  kl  of  GaSb  [black 
circles  (Ref.  56)],  GaAs  [green  circles  (Ref.  47)  and  green  triangles 
(Ref.  48)],  and  GaP  [red  circles  (Ref.  5 1)  and  red  triangles  (Ref.  52)]. 

B.  Gallium-Y  compounds 

We  previously  presented  results  for  K\  and  P  for  wurtzite 
GaN  and  compared  those  results  with  the  calculated  P  for 
GaAs,  GaSb,  GaP,  and  cubic  GaN.16  Here  we  give  k\  for 
these  cubic  compounds.  Figure  5  shows  the  calculated  /cnaturai 
versus  temperature  for  GaSb  (black  dashed  curve),  GaAs 
(green  dashed  curve),  GaP  (red  dashed  curve),  and  cubic 
GaN  (blue  dashed  curve).  Also  shown  are  /cpure  for  GaP  (red 
solid  curve)  and  cubic  GaN  (blue  solid  curve).  Experimental 
data  is  shown  for  /cl  of  GaSb  (black  circles56),  GaAs  (green 
circles47  and  green  triangles48),  and  GaP  (red  circles51  and  red 
triangles52).  For  GaN  the  lattice  constant  determined  by  energy 
minimization  was  increased  ~1%  to  give  better  agreement 
with  the  experimental  phonon  dispersion  as  discussed  in 
Ref.  16  and  discussed  below  for  the  indium- V  materials.  We 
used  the  lattice  constants  determined  by  energy  minimization 
for  GaAs,  GaP,  and  GaSb,  which  gave  good  agreement  for  the 
phonon  dispersions  (see  Figs.  14-16  in  the  Appendix). 

GaN  has  the  largest  frequency  gap  between  acoustic  and 
optic  phonons  and  the  largest  acoustic  frequency  scale  of 
the  gallium-V  materials  (see  Figs.  14-17  in  the  Appendix). 
These  properties  give  higher  phonon  velocities  and  weaker 
anharmonic  phonon  scattering  in  GaN  and  thus  significantly 
larger  k\  than  the  other  materials  in  Fig.  5.  /cpurc  for  GaN 
is  over  two  times  larger  than  /cpure  for  GaP  over  the  entire 
temperature  range.  We  also  note  that  the  full  BTE  solution 
for  /cpure  of  GaN  is  13%  higher  than  the  RTA  result  at  T  = 
300  K.  For  the  other  materials  the  difference  is  around  5%  or 
less  (see  Table  I).  A  disparity  exists  between  the  calculated 
^natural  for  GaP  and  the  experimental  results  from  Ref.  52, 
the  calculated  result  being  20%  higher  at  room  temperature. 
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temperature  (K) 

FIG.  6.  (Color  online)  Calculated  isotope  effect  P  vs  temperature 
for  cubic  GaN  (solid  black  curve),  wurtzite  GaN  ( dashed  black  curve), 
GaP  (green  curve),  GaSb  (orange  curve),  diamond  (gold  curve),  Ge 
(red  curve),  AlSb  (brown  curve),  and  3C-SiC  (blue  curve). 

The  samples  studied  in  Ref.  52  were  believed  to  have  had 
appreciable  amounts  of  impurities  which  may  be  the  cause  of 
the  differences.  The  calculated  /cnaturai  for  GaAs  and  GaSb  are 
in  good  agreement  with  experiment. 

Due  to  the  large  frequency  gap  and  frequency  scale  in 
GaN,  the  anharmonic  phonon  scattering  is  relatively  weak 
and  phonon-isotope  scattering  is  relatively  more  important  for 
limiting  Kt  than  in  materials  with  stronger  phonon-phonon 
scattering.  The  isotope  effect  in  GaN  is  very  large,  with  P  — 
68%  at  T  —  300  K  and  increases  with  decreasing  temperature 
as  the  anharmonic  phonon  scattering  gets  weaker.  Figure  6 
shows  previously  published16  calculated  isotope  effects  P 
versus  temperature  for  cubic  GaN  (solid  black  curve),  wurtzite 
GaN  (dashed  black  curve),  GaP  (green  curve),  and  GaSb 
(orange  curve)  with  the  addition  of  the  calculated  P  for 
diamond11  (gold  curve),  Ge  (red  curve),  AlSb  (brown  curve), 
and  3C-SiC  (blue  curve).  The  room  temperature  P  for  all 
the  materials  studied  in  this  work  are  also  given  in  Table  I. 
At  T  =  300  K,  GaN  has  the  largest  isotope  effect  which 
can  be  understood  in  terms  of  the  large  Ga  isotope  mixture, 
the  large  frequency  gap,  and  the  high  frequency  scale  in 
GaN.16  Though  diamond  has  a  smaller  isotopic  mass  variance 
than  GaN  and  no  frequency  gap  between  acoustic  and  optic 
phonons,  it  also  has  a  large  isotope  effect,  P  =  51%  at  T  — 
300  K.  Anharmonic  phonon  coupling  is  very  weak  in  diamond 
at  T  =  300  K,  making  phonon-isotope  scattering  important. 
The  weak  anharmonic  coupling  arises  from  a  small  unit  cell 
mass  and  strong  covalent  bonding  which  gives  diamond  a  very 
high  acoustic  frequency  scale,  ~3.5  times  higher  than  GaN 
(compare  Fig.  17  in  the  appendix  of  this  work  and  Fig.  1  of 
Ref.  9).  These  properties  give  diamond  a  high  K\  and  P.9A  1 

Despite  having  a  smaller  isotopic  mass  variance  and 
relatively  small  frequency  scale  (roughly  half  that  of  GaN), 
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AlSb  has  P  —  36%  at  T  =  300  K.  The  large  isotope  effect  and 
its  temperature  dependence  are  a  consequence  of  significant 
optic  phonon  contributions  to  kl  discussed  above.  The  weak 
anharmonic  scattering  leads  to  significantly  increased  optic 
phonon  lifetimes  and  larger  contributions  to  /rpuie  compared  to 
other  materials.  The  isotope  scattering  rates  scale  roughly  as 
or  Dim), 26  where  l)(o>)  is  the  phonon  density  of  states  which 
is  large  for  the  relatively  flat  optic  branches.  Thus,  isotope 
scattering  strongly  suppresses  the  optic  phonon  contributions 
to  /Cnaturah  leading  to  a  large  percentage  isotope  effect  P. 
Unlike  the  other  curves  in  Fig.  6,  the  P  for  AlSb  has  a  peak  at 
175  K  and  decreases  with  decreasing  temperature  below  the 
peak.  As  the  temperature  decreases  below  the  peak,  the  optic 
phonons  become  less  thermally  populated  and  contribute 
less  to  kl  than  the  acoustic  phonons  despite  weakening 
anharmonic  scattering.  Since  the  acoustic  phonons  have  lower 
frequencies  and  are  less  affected  by  the  isotope  scattering 
than  the  optic  phonons,  P  is  reduced. 

C.  Indium-V  compounds 

As  discussed  in  Sec.  Ill,  the  LDA  approach  to  calculate 
the  harmonic  and  anharmonic  IFCs  tends  to  overbind  atoms,57 
and  corrections  for  zero-point  and  nonzero  temperature  atomic 
motion58  account  for  only  a  small  fraction  of  the  difference.  For 
most  of  the  systems  examined  in  this  work,  excellent  agree¬ 
ment  with  the  experimental  phonon  dispersions  (see  figures 
in  the  Appendix)  is  achieved  in  spite  of  this  underestimation 
of  the  lattice  constants.  However,  similar  to  previous  work 
calculating  phonon  dispersions  and  k\  of  GaN,  we  find  the 
calculated  optic  phonon  frequencies  for  InP,  InAs,  and  InSb  are 
higher  than  the  experimental  dispersions,  leading  to  artificially 
weakened  acoustic-optic  phonon  scattering. 

Energy  minimization  gives  lattice  constants  <7inp  =  5.79  A, 
aInAs  =  5.97  A,  and  «lnSb  =  6.39  A  while  measured  values  are 
1.43%,  1.53%,  and  1.40%  higher,  respectively.  Further,  the 
optic  phonon  branches  given  by  these  lattice  constants  (red 
curves  in  Figs.  18-20  in  the  Appendix)  lie  above  the  experi¬ 
mental  phonon  dispersions.  The  black  curves  in  these  figures 
correspond  to  the  calculated  phonon  dispersions  using  lattice 
constants  that  were  increased  by  1.5%  for  InP  and  1.0%  for 
InAs  and  InSb.  Similar  to  GaN,  the  increased  lattice  constants 
for  the  indium-V  systems  lower  the  high-lying  optic  phonon 
branches  and  give  better  agreement  with  the  experimental  data. 
The  calculated  Arnaturai  given  by  the  different  lattice  constants 
for  each  indium-V  compound  are  compared  in  Fig.  7. 

Figure  7  shows  the  calculated  a:l  versus  temperature  for 
the  indium-V  compounds  with  experimental  data  for  InP 
(black  circles,59  squares,60  diamond,61  and  triangle62),  InAs 
(red  circles,63  squares,64  diamond,65  and  triangles66),  and 
InSb  (blue  circles,56  squares,6  and  triangles68).  The  isotope 
effect  is  small  in  the  indium-V  compounds  so  we  show 
only  the  calculated  ATnaturai  in  Fig-  7.  The  dashed  black  curve 
corresponds  to  /rnaturai  for  InP  with  harmonic  and  anharmonic 
IFCs  determined  using  the  lattice  constant  from  minimization 
of  energy,  alnp  =  5.79  A.  The  solid  black  curve  corresponds 
to  /cnaturai  for  InP  using  the  corrected  lattice  constant  ainp  = 
5.87  A.  The  calculated  /cnaturai  of  InAs  is  given  by  the  dashed 
[aInAs  =  5.97  A  (energy  minimization)]  and  solid  [c/inAs  = 
6.03  A  (corrected)]  red  curves  and  /cnaturai  of  InSb  is  given 


temperature  (K) 


FIG.  7.  (Color  online)  Calculated  k l  vs  temperature  for  the 
indium-V  compounds  with  experimental  data  for  InP  [black  circles 
(Ref.  59),  squares  (Ref.  60),  diamond  (Ref.  61),  and  triangle 
(Ref.  62)],  InAs  [red  circles  (Ref.  63),  squares  (Ref.  64),  diamond 
(Ref.  65),  and  triangles  (Ref.  66)],  and  InSb  [blue  circles  (Ref.  56), 
squares  (Ref.  67),  and  triangles  (Ref.  68)].  The  dashed  black  curve 
corresponds  to  /c„aturai  for  InP  with  harmonic  and  anharmonic  IFCs 
determined  using  the  lattice  constant  from  energy  minimization, 
a[nP  =  5.79  A.  The  solid  black  curve  corresponds  to  /cna tural  using 
the  corrected  lattice  constant,  a\np  =  5.87  A.  Calculated  at„ aturai  of 
InAs  is  given  by  the  dashed  [flinAs  —  5.97  A  (energy  minimization)] 
and  solid  [a^s  =  6.03  A  (corrected)]  red  curves  and  /cnaturai  of  InSb 
is  given  by  dashed  [flinsb  =  6.39  A  (energy  minimization)]  and  solid 
[flinSb  =  6.45  A  (corrected)]  blue  curves. 

by  dashed  [«insb  =  6.39  A  (energy  minimization)]  and  solid 
[nInSb  =  6.45  A  (corrected)]  blue  curves. 

For  all  indium-V  compounds  the  /^natural  given  by  the 
corrected  lattice  constants  (solid  curves)  lie  below  the  Arnaturai 
given  by  the  lattice  constants  determined  by  energy  mini¬ 
mization  (dashed  curves).  With  the  increased  lattice  constants 
the  higher-lying  optic  modes  are  shifted  to  lower  frequencies 
which  provide  stronger  scattering  channels  for  the  heat¬ 
carrying  acoustic  phonons  and  thus  reduced  atl-  At  T  —  300 
K,  /^natural  determined  by  the  corrected  lattice  constants  are 
3%,  8%,  and  5%  lower  in  InP,  InAs,  and  InSb,  respectively, 
which  gives  somewhat  better  agreement  with  experiment.  In 
general,  the  calculated  Arnaturai  for  InSb  is  in  good  agreement 
with  experiment.  The  calculated  /cnatura i  for  InAs  is  higher  than 
experiment  throughout  the  temperature  range,  suggesting  that 
a  crystal  with  fewer  defects  could  lead  to  modest  improvements 
in  atl-  /^natural  for  InP  is  in  agreement  with  the  experimental  data 
point  at  T  —  300  K  given  in  Ref.  62,  but  is  higher  than  the  other 
experimental  data  throughout  the  temperature  range.  We  also 
note  that  the  RTA  solution  to  the  BTE  is  a  good  approximation 
in  the  indium-V  compounds  and  gives  arta  only  a  few  percent 
below  the  full  BTE  solution. 

atl  in  the  compound  semiconductors  depends  critically  on 
the  anharmonic  three-phonon  scattering  rates  determined  from 
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AB  INITIO  THERMAL  TRANSPORT  IN  COMPOUND  . . . 


perturbation  theory: 
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No  is  the  number  of  unit  cells  in  the  crystal  and  the  ± 
signs  correspond  to  the  two  types  of  possible  three-phonon 
processes  that  satisfy  conservation  of  crystal  momentum  and 
conservation  of  energy.1  The  top  term  in  the  brackets  is  for 
the  +  process,  and  the  bottom  term  for  the  —  process.  The 
scattering  matrix  elements  depend  on  the  lowest 

order  anharmonic  IFCs,  Q>apY(lic,l'K',  1"k").  We  note  that  the 
three-phonon  scattering  rates  depend  on  both  harmonic  and  an¬ 
harmonic  properties  of  each  system.  Stiff  interatomic  bonding 
in  systems  such  as  diamond,  SiC,  AIN,  and  GaN  not  only  give 
large  acoustic  phonon  velocities  and  high  frequency  scales,  but 
also  larger  anharmonic  IFCs.  Of  all  the  systems  we  have  con¬ 
sidered,  diamond  has  the  largest  kl  and  the  largest  anharmonic 
IFCs,  which  tend  to  increase  the  three-phonon  scattering  rates. 
However,  the  higher  overall  phonon  frequencies  of  diamond 
enter  the  denominator  of  Eq.  (5)  and  counteract  the  effects 
of  the  larger  anharmonic  IFCs  to  give  lower  scattering  rates. 
Conversely,  InSb  has  lower  overall  phonon  frequencies,  has 
generally  smaller  anharmonic  IFCs,  and  has  the  lowest  of 
the  compound  systems.  Further,  the  intrinsic  phonon  transport 
lifetimes  also  depend  on  the  phase  space  for  allowed  scattering 
processes  given  by  the  conservation  conditions  and  the 
coupling  of  nonequilibrium  phonon  modes  through  the  BTE. 


V.  SILICON  CARBIDE  AND  ALUMINUM  NITRIDE 

Here  we  apply  the  first-principles  BTE  approach  to 
calculating  k\  for  technologically  important  semiconductors, 
cubic  3C-SiC  and  wurtzite  AIN.  To  determine  the  anharmonic 
IFCs  for  AIN  we  used  a  108  atom  supercell  with  an  interaction 
cutoff  to  include  fourth  nearest  neighbors  similar  to  that 
used  previously  for  GaN.16  For  cubic  SiC  we  used  216  atom 
supercells  with  an  interaction  cutoff  to  include  third  nearest 
neighbors,  similar  to  the  cubic  compounds  discussed  above. 
For  each  system  we  enforced  translational  invariance  via  the 
X 1  minimization  procedure.  Figure  8  shows  the  calculated 
/Cpure  (solid  black  curve)  and  /rnaturai  (dashed  black  curve)  for 
3C-SiC  compared  with  experimental  data  given  by  circles.69 
The  calculated  in-plane  and  out-of-plane  (c-axis)  /cpure  for 
AIN  are  given  by  the  red  and  green  curves,  respectively.  Open 
green  circles  are  the  experimental  atl  along  the  c  axis  for 
AIN. 70  The  isotope  mixtures  for  A1  and  N  are  negligible  so 
we  do  not  give  /cnaturai  for  AIN. 

Light  atomic  masses  and  stiff  atomic  bonds  for  the  Si 
and  C  atoms  give  SiC  large  acoustic  phonon  velocities  and 
a  high  phonon  frequency  scale,  though  significantly  smaller 
than  those  of  diamond.  The  relatively  large  Si  to  C  mass 
ratio  produces  an  appreciable  frequency  gap  between  acoustic 
and  optic  branches  as  seen  in  Fig.  21  in  the  Appendix. 
These  properties  make  SiC  a  high  atl  material  with  calculated 
Natural  =  480  Wm”1  KT1  and  ^pure  =  570  Wm"1  K”1  and  a 
percentage  isotope  effect  of  P  —  20%  at  300  K.  We  note  that 
A:pUre  of  cubic  SiC  is  still  six  times  smaller  than  that  of  diamond. 
AIN  is  also  a  relatively  high  k l  material  with  calculated  Apurt  = 


temperature  (K) 


FIG.  8.  (Color  online)  Calculated  Arpure  (solid  black  curve)  and 
^natural  (dashed  black  curve)  for  3C-SiC  with  experimental  data  given 
by  circles  (Ref.  69).  The  calculated  in-plane  and  out-of-plane  (c-axis) 
A:pUre  for  AIN  are  given  by  the  red  and  green  curves,  respectively,  with 
experimental  data  for  c-axis  kl  (Ref.  70). 

322  W  m_1  K  1  at  T  —  300  K.  For  each  system  the  calculated 
kl  gives  good  agreement  with  the  experimental  data  around 
and  above  room  temperature.  At  lower  temperatures  where 
scattering  from  defects  and  boundaries  becomes  more  impor¬ 
tant,  each  calculated  a:l  is  predictably  higher  than  experiment. 

We  found  previously  that  the  wurtzite  and  cubic  GaN 
structures  have  k\  that  are  similar  and  larger  than  the  thermal 
conductivities  of  the  cubic  GnX  compounds  (X  =  P,  As, 
and  Sb).16  Similarly,  wurtzite  AIN  has  significantly  higher 
a:l  than  the  cubic  A1 A  compounds  discussed  in  the  previous 
section  (see  Table  I).  AIN  has  the  largest  acoustic  phonon 
frequency  scale  and  the  highest  acoustic  velocities  of  all  the 
A1X  materials.  Thus,  at  100  K,  the  in-plane  A:pure  of  AIN  is  the 
largest  of  these  materials,  over  five  times  larger  than  that  of  A1P. 
AIN  does  not  have  an  appreciable  gap  between  the  acoustic 
and  optic  branches  (see  Fig.  22  in  the  Appendix),  and  thus  with 
increasing  temperature  the  kl  of  AIN  decreases  rapidly  as  aao 
scattering  becomes  stronger.  At  300  K,  Arpure  is  still  more  than 
2.5  times  larger  than  that  for  the  cubic  A1X  compounds.  The 
decay  of  optic  phonons  into  the  heat-carrying  acoustic  phonons 
was  previously  found  to  increase  a:l  in  AIN.54  Here,  we  show 
that  the  coupling  of  acoustic  and  optic  phonons  reduces  the  a:l 
for  the  systems  studied  here.  The  coupling  of  different  phonon 
modes  through  the  full  solution  to  the  BTE  plays  a  significant 
role  in  determining  atl  of  AIN.  The  single  mode  relaxation 
time  approximation  gives  atrta  of  AIN  14%  smaller  than  that 
given  by  the  full  BTE  solution. 


VI.  SUMMARY  AND  CONCLUSIONS 

We  have  presented  an  accurate  and  predictive  first- 
principles  method  for  calculating  the  lattice  thermal  con¬ 
ductivity  atl  for  a  range  of  semiconducting  systems.  This 
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method  uses  harmonic  and  anharmonic  interatomic  forces 
from  density  functional  theory  calculations  as  inputs  to  the 
linearized  phonon  Boltzmann  transport  equation  for  which  an 
exact  numerical  solution  is  obtained.  The  calculations  of  the 
anharmonic  force  constants  and  their  role  in  determining  the 
thermal  conductivities  of  a  variety  of  systems  was  discussed. 
This  first-principles  approach  was  tested  on  well-studied 
systems  (Si,  Ge,  and  GaAs)  and  excellent  agreement  with 
measured  k\  data  was  found.  We  used  this  approach  to  examine 
kl  in  aluminum- V,  gallium- V,  and  indium-V  compounds  as 
well  as  the  technologically  important  materials,  SiC  and  AIN. 
The  upper  bound  to  the  intrinsic  K{  as  limited  by  anharmonic 
phonon  scattering  was  calculated  for  each  system,  and  the  role 
of  isotope  scattering  was  examined.  The  interplay  between 
harmonic  and  anharmonic  crystal  properties  in  determining  the 
intrinsic  phonon-phonon  scattering  rates  was  also  discussed. 

In  general,  the  calculated  Arnatumi  is  in  good  agreement  with 
the  available  experimental  data  for  most  materials.  For  some 
compounds  such  as  GaP,  AlSb,  InAs,  and  InP  the  calculated 
^natural  lies  above  the  measured  data,  suggesting  that  improved 
crystal  growth  to  reduce  point  defects  and  other  extrinsic 
scattering  mechanisms  could  lead  to  enhancements  in  /cL.  We 
have  found  high  intrinsic  k l  values  for  3C-SiC,  AIN,  and  cubic 
GaN  ranging  from  around  300  to  600  W  m  1  K'1  at  T  —  300 
K,  while  cubic  GaN  and  AlSb  showed  large  isotope  effects.  The 
frequency  gap  between  acoustic  and  optic  phonons  plays  an 
important  role  in  reducing  the  acoustic-optic  phonon  scattering 
and  in  determining  k\  for  compounds  with  large  mass  differ¬ 
ences  between  the  cations  and  the  anions.  A  large  acoustic- 
optic  frequency  gap  in  AlSb  leads  to  significant  optic  phonon 
contributions  to  K\  and  to  a  large  calculated  isotope  effect. 

The  microscopic,  first-principles  method  presented  in  this 
work  is  a  powerful  tool  for  calculating  and  examining  kl  which 
is  predictive  and  transferable  to  a  wide  range  of  systems. 


package.34,35  For  most  cases  we  used  norm-conserving  pseu¬ 
dopotentials  with  80  Ryd  plane-wave  cutoffs  and  6x6x6 
Monkhorst-Pack  A; -point  meshes  for  electronic  and  phonon 
calculations.  All  calculations  were  made  with  pseudopotentials 
from  the  QUANTUM  ESPRESSO  website'4  and  the  names  are 
given  in  the  figure  captions. 

Here  we  give  the  phonon  dispersions  in  the  high-symmetry 
directions  for  all  of  the  materials  considered  in  this  work. 
Available  experimental  data  is  included  for  comparison  with 
the  calculated  curves  (see  Figs.  9-22). 


FIG.  9.  Calculated  phonon  dispersion  for  Si  in  the  indicated  high 
symmetry  directions  (black  curves).  The  pseudopotential  used  in  the 
LDA/DFPT  calculation  was  Si.pz-vbc.UPF  (Ref.  34).  Experimental 
data  are  given  by  black  circles  (Ref.  71). 
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APPENDIX 

Given  a  set  of  harmonic  IFCs,  <t>ap(lK,l'i<:'),  the  phonon  fre¬ 
quencies  are  determined  by  diagonalization  of  the  dynamical 
matrix 

DKaUq)  =  1  _  V  4\tf(0 k/k'V5  *"  ,  (Al) 

for  a  given  wave  vector  q.  Here,  mK  is  the  isotope  averaged 
mass  of  the  /cth  atom,  R/  is  the  lattice  vector  locating  the  / th 
unit  cell,  and  a  and  /3  are  Cartesian  directions.  The  harmonic 
IFCs  for  each  system  were  calculated  within  the  framework 
of  DFPT  and  LDA  using  the  plane-wave  QUANTUM  ESPRESSO 
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FIG.  10.  Calculated  phonon  dispersion  for  Ge  in  the  high  sym¬ 
metry  directions  (black  curves).  The  pseudopotential  used  in  the 
LDA/DFPT  calculation  was  Ge.pz-bhs.UPF  (Ref.  34).  Experimental 
data  are  given  by  black  circles  (Ref.  71). 
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FIG.  11.  Calculated  phonon  dispersion  for  A1P  in  the  high 
symmetry  directions  (black  curves).  The  pseudopotentials  used  in 
the  LDA/DFPT  calculation  were  Al.pz-vbc.UPF  and  P.pz-bhs.UPF 
(Ref.  34).  Experimental  data  are  given  by  black  circles  (Ref.  72). 


FIG.  13.  Calculated  phonon  dispersion  for  AlSb  in  the  high 
symmetry  directions  (black  curves).  The  pseudopotentials  used  in 
the  LDA/DFPT  calculation  were  Al.pz-vbc.UPF  and  Sb.pz-bhs.UPF 
(Ref.  34).  Experimental  data  are  given  by  black  circles  (Ref.  74). 


FIG.  12.  Calculated  phonon  dispersion  for  AlAs  in  the  high 
symmetry  directions  (black  curves).  The  pseudopotentials  used  in 
the  LDA/DFPT  calculation  were  Al.pz-vbc.UPF  and  As.pz-bhs.UPF 
(Ref.  34).  Experimental  data  are  given  by  black  circles  (Ref.  73). 


FIG.  14.  Calculated  phonon  dispersion  for  GaP  in  the  high 
symmetry  directions  (black  curves).  The  pseudopotentials  used  in 
the  LDA/DFPT  calculation  were  Ga.pz-bhs.UPF  and  P.pz-bhs.UPF 
(Ref.  34).  Experimental  data  are  given  by  black  circles  (Ref.  75). 


165201-11 


L.  LINDSAY,  D.  A.  BROIDO,  AND  T.  L.  REINECKE 


PHYSICAL  REVIEW  B  87,  165201  (2013) 


FIG.  15.  Calculated  phonon  dispersion  for  GaAs  in  the  high 
symmetry  directions  (black  curves).  The  pseudopotentials  used  in 
the  LDA/DFPT  calculation  were  Ga.pz-bhs.UPF  and  As.pz-bhs.UPF 
(Ref.  34).  Experimental  data  are  given  by  black  circles  (Ref.  76). 


FIG.  17.  Calculated  phonon  dispersion  for  cubic  GaN  in  the  high 
symmetry  directions  (black  curves).  The  pseudopotentials  used  in 
the  LDA/DFPT  calculation  were  Ga.pz-bhs.UPF  and  N.pz-vbc.UPF 
(Ref.  34).  We  note  that  the  calculated  dispersion  was  obtained  using 
a  1%  increase  to  the  lattice  constant  from  energy  minimization  as 
discussed  in  Ref.  16. 


scaled  wave  vector 


FIG.  16.  Calculated  phonon  dispersion  for  GaSb  in  the  high 
symmetry  directions  (black  curves).  The  pseudopotentials  used  in 
the  LDA/DFPT  calculation  were  Ga.pz-bhs.UPF  and  Sb.pz-bhs.UPF 
(Ref.  34).  Experimental  data  are  given  by  black  circles  (Ref.  77). 


FIG.  18.  (Color  online)  Calculated  phonon  dispersion  for  InP 
in  the  high  symmetry  directions.  Red  curves  give  dispersion  using 
a  lattice  constant  determined  by  energy  minimization  and  black 
curves  give  dispersion  using  this  lattice  constant  with  a  1.5% 
increase.  The  pseudopotentials  used  in  the  LDA/DFPT  calculation 
were  In.pz-n-bhs.UPF  and  P.pz-bhs.UPF  (Ref.  34).  We  note  that  the 
indium  pseudopotential  included  a  core  correction.  Experimental  data 
are  given  by  black  circles  (Ref.  78). 
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FIG.  19.  (Color  online)  Calculated  phonon  dispersion  for  InAs 
in  the  high  symmetry  directions.  Red  curves  give  dispersion  using  a 
lattice  constant  determined  by  energy  minimization  and  black  curves 
give  dispersion  using  this  lattice  constant  with  a  1%  increase.  The 
pseudopotentials  used  in  the  LDA/DFPT  calculation  were  In.pz-n- 
bhs.UPF  and  As.pz-bhs.UPF  (Ref.  34).  We  note  that  the  indium 
pseudopotential  included  a  core  correction.  Experimental  data  are 
given  by  black  circles  (Ref.  79)  and  black  squares  (Ref.  80). 
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FIG.  20.  (Color  online)  Calculated  phonon  dispersion  for  InSb 
in  the  high  symmetry  directions.  Red  curves  give  dispersion  using  a 
lattice  constant  determined  by  energy  minimization  and  black  curves 
give  dispersion  using  this  lattice  constant  with  a  1%  increase.  The 
pseudopotentials  used  in  the  LDA/DFPT  calculation  were  In.pz-n- 
bhs.UPF  and  Sb.pz-bhs.UPF  (Ref.  34).  We  note  that  the  indium 
pseudopotential  included  a  core  correction.  Experimental  data  are 
given  by  black  circles  (Ref.  81). 
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FIG.  21.  Calculated  phonon  dispersion  for  cubic  3C-SiC  in  the 
high  symmetry  directions  (black  curves).  The  pseudopotentials  used 
in  the  LDA/DFPT  calculation  were  Si.pz-vbc.UPF  and  C.pz-vbc.UPF 
(Ref.  34).  Experimental  data  are  given  by  black  circles  (Ref.  82). 


FIG.  22.  Calculated  phonon  dispersion  for  wurtzite  AIN  in  the 
high  symmetry  directions  (black  curves).  The  pseudopotentials 
used  in  the  LDA/DFPT  calculation  were  Al.pz-vbc.UPF  and  N.pz- 
vbc.UPF  (Ref.  34).  Experimental  data  are  given  by  black  circles 
(Ref.  83). 


165201-13 


L.  LINDSAY,  D.  A.  BROIDO,  AND  T.  L.  REINECKE 


PHYSICAL  REVIEW  B  87,  165201  (2013) 


1 J.  M.  Ziman,  Electrons  and  Phonons  (Oxford  University  Press, 
London,  1960). 

2P.  G.  Klemens,  Solid  State  Phys.  7,  1  (1958). 

3J.  Callaway,  Phys.  Rev.  113,  1046  (1958). 

4M.  G.  Holland,  Phys.  Rev.  132,  2461  (1963). 

5P.  G.  Klemens,  J.  Wide  Bandgap  Mater.  7,  332  (2000). 

6R.  E.  Peierls,  Quantum  Theory  of  Solids  (Oxford  University  Press, 
London,  1955). 

7D.  A.  Broido,  M.  Malorny,  G.  Birner,  N.  Mingo,  and  D.  A.  Stewart, 
Appl.  Phys.  Lett.  91,  231922  (2007). 

8J.  Garg,  N.  Bonini,  B.  Kozinsky,  and  N.  Marzari,  Phys.  Rev.  Lett. 
106,045901  (2011). 

9A.  Ward,  D.  A.  Broido,  D.  A.  Stewart,  and  G.  Deinzer,  Phys.  Rev. 
B  80,  125203  (2009). 

10Wu  Li,  N.  Mingo,  L.  Lindsay,  D.  A.  Broido,  D.  A.  Stewart,  and  N. 

A.  Katcho,  Phys.  Rev.  B  85,  195436  (2012). 
nD.  A.  Broido,  L.  Lindsay,  and  A.  Ward,  Phys.  Rev.  B  86,  115203 
(2012). 

12T.  Luo,  J.  Garg,  J.  Shiomi,  K.  Esfarjani,  and  G.  Chen,  Europhys. 
Lett.  101,  16001  (2013). 

13X.  Tang  and  J.  Dong,  Proc.  Natl.  Acad.  Sci.  USA  107, 4539  (2010). 
14J.  Shiomi,  K.  Esfarjani,  and  G.  Chen,  Phys.  Rev.  B  84,  104302 
(2011). 

15Z.  Tian,  J.  Garg,  K.  Esfarjani,  T.  Shiga,  J.  Shiomi,  and  G.  Chen, 
Phys.  Rev.  B  85,  184303  (2012). 

16L.  Lindsay,  D.  A.  Broido,  and  T.  L.  Reinecke,  Phys.  Rev.  Lett.  109, 
095901  (2012). 

l7Wu  Li,  L.  Lindsay,  D.  A.  Broido,  D.  A.  Stewart,  and  N.  Mingo, 
Phys.  Rev.  B  86,  174307  (2012). 

18P.  Hohenberg  and  W.  Kohn,  Phys.  Rev.  136,  B864  (1964). 

19W.  Kohn  and  L.  J.  Sham,  Phys.  Rev.  140,  A1133  (1965). 

20S.  Baroni,  S.  Gironcoli,  A.  D.  Corso,  and  P.  Giannozzi,  Rev.  Mod. 
Phys.  73,515  (2001). 

21 K.  Esfarjani,  G.  Chen,  and  H.  T.  Stokes,  Phys.  Rev.  B  84,  085204 

(2011). 

22L.  Lindsay,  D.  A.  Broido,  andN.  Mingo,  Phys.  Rev.  B  82, 161402(R) 

(2010). 

23L.  Lindsay,  D.  A.  Broido,  and  N.  Mingo,  Phys.  Rev.  B  83,  235428 

(2011). 

24L.  Lindsay  and  D.  A.  Broido,  Phys.  Rev.  B  85,  035436  (2012). 

25  We  note  that  /cl  of  wurtzite  GaN  and  AIN  have  small  in-plane 
anisotropies  of  around  a  few  percent  around  room  temperature  and 
higher.  We  define  tcm  =  (kxx  +  Kyy)/ 2,  where  the  hexagonal  plane 
lies  in  the  x-y  plane  and  the  c  axis  is  along  the  z  direction. 

26S.  I.  Tamura,  Phys.  Rev.  B  30,  849  (1984). 

27  A  list  of  the  mass  variance  parameters  for  naturally  occurring 
isotope  abundances  of  the  materials  studied  in  this  work  is  found  in 
Table  I. 

28M.  Omini  and  A.  Sparavigna,  Phys.  Rev.  B  53,  9064  (1996). 

29M.  Omini  and  A.  Sparavigna,  Nuovo  Cimento  Soc.  Ital.  Fis..  D  19, 
1537  (1997). 

30A.  Ward  and  D.  A.  Broido,  Phys.  Rev.  B  81,  085205  (2010). 

31  J.  H.  Seol,  I.  Jo,  A.  L.  Moore,  L.  Lindsay,  Z.  H.  Aitken,  M.  T.  Pettes, 
X.  Li,  Z.  Yao,  R.  Huang,  D.  A.  Broido,  N.  Mingo,  R.  S.  Ruoff,  and 
L.  Shi,  Science  328,  213  (2010). 

32L.  Lindsay,  D.  A.  Broido,  and  N.  Mingo,  Phys.  Rev.  B  82,  1 15427 

(2010). 

33L.  Lindsay,  D.  A.  Broido,  and  N.  Mingo,  Phys.  Rev.  B  80,  125407 
(2009). 


34S.  Baroni  et  al.,  http://www.quantum-espresso.org. 

35P.  Gianozzi  et  al.,  J.  Phys.:  Condens.  Matter  21,  395502  (2009). 

36For  the  wurtzite  structures  the  internal  parameter  u,  which  defines 
the  relative  position  of  the  two  hexagonal  close-packed  sublattices, 
was  chosen  so  that  the  Hellmann-Feynman  forces  on  the  atoms 
were  negligible. 

37K.  Esfarjani  and  H.  T.  Stokes,  Phys.  Rev.  B  77,  1441 12  (2008). 

38G.  Leibfried  and  W.  Ludwig,  Solid  State  Phys.  12,  275  (1961). 

39G.  P.  Srivastava,  The  Physics  of  Phonons  (Taylor  and  Francis,  New 
York,  1990). 

40W.  H.  Press,  S.  A.  Teukolsky,  W.  T.  Vetterling,  and  B.  P. 
Flannery,  Numerical  Recipes  in  Fortran  (Cambridge  University 
Press,  Cambridge,  UK,  1992). 

41N.  Mingo,  D.  A.  Stewart,  D.  A.  Broido,  and  D.  Srivastava,  Phys. 
Rev.  B  77,  033418  (2008). 

42A.  V.  Inyushkin,  A.  N.  Taldenkov,  A.  M.  Gibin,  A.  V.  Gusev,  and 
H.-J.  Pohl,  Phys.  Status  Solidi  C  1,  2995  (2004). 

43T.  Ruf,  R.  W.  Henn,  M.  Asen-Palmer,  E.  Gmelin,  M.  Cardona,  H.-J. 
Pohl,  G.  G.  Devyatych,  and  P.  G.  Sennikov,  Solid  State  Commun. 
115,  243  (2000). 

44We  note  that  of  isotopically  enriched  Si  was  also  reported  in 
Ref.  43  (not  shown  here),  however,  the  room  temperature  values 
were  different  from  measurements  reported  in  Ref.  45  on  the  same 
sample:  T.  Ruf,  R.  W.  Henn,  M.  Asen-Palmer,  E.  Gmelin,  M. 
Cardona,  H.-J.  Pohl,  G.  G.  Devyatych,  and  P.  G.  Sennikov,  Solid 
State  Commun.  127,  257  (2003). 

4iR.  K.  Kremer,  K.  Graf,  M.  Cardona,  G.  G.  Devyatykh,  A.  V.  Gusev, 
A.  M.  Gibin,  A.  V.  Inyushkin,  A.  N.  Taldenkov,  and  H.-J.  Pohl, 
Solid  State  Commun.  131,  499  (2004). 

46V.  I.  Ozhogin,  A.  V.  Inyushkin,  A.  N.  Taldenkov,  A.  V.  Tikhomirov, 
G.  E.  Popov,  E.  Haller,  and  K.  Itoh,  JETP  Lett.  63,  490  (1996). 

47R.  O.  Carlson,  G.  A.  Slack,  and  S.  J.  Silverman,  J.  Appl.  Phys.  36, 
505  (1965). 

48 A.  V.  Inyushkin,  A.  N.  Taldenkov,  A.  Y.  Yakubovsky,  A.  V.  Markov, 
L.  Moreno-Garsia,  and  B.  N.  Sharonov,  Semicond.  Sci.  Technol. 
18,  685  (2003). 

49G.  Deinzer,  G.  Birner,  and  D.  Strauch,  Phys.  Rev.  B  67.  144304 
(2003). 

50G.  Deinzer,  M.  Schmitt,  A.  P.  Mayer,  and  D.  Strauch,  Phys.  Rev.  B 
69,  014304  (2004). 

51E.  F.  Steigmeier  and  I.  Kudman,  Phys.  Rev.  141,  767  (1966). 

52V.  M.  Muzhdaba,  A.  Y.  Nashel’skii,  P.  V.  Tamarin,  and  S.  S.  Shalyt, 
Sov.  Phys.  Solid  State  10,  2265  (1969). 

53M.  A.  Afromowitz,  J.  Appl.  Phys.  44,  1292  (1973). 

54M.  Kazan,  S.  Pereira,  and  M.  R.  Correia,  Phys.  Rev.  B  77, 1 80302(R) 
(2008). 

55L.  Lindsay  and  D.  A.  Broido,  J.  Phys.:  Condens.  Matter  20,  165209 
(2008). 

56M.  G.  Holland,  Phys.  Rev.  134,  A471  (1964). 

57 P.  Haas,  F.  Tran,  and  P.  Blaha,  Phys.  Rev.  B  79,  085104  (2009). 

58P.  B.  Allen,  Philos.  Mag.  B  70,  527  (1994). 

59S.  A.  Aliev,  A.  Ya.  Nashel'skii,  and  S.  S.  Shalyt,  Sov.  Phys.  Solid 
State  7,  1287  (1965). 

60S.  Adachi,  Properties  and  Applications  of  Indium  Phosphide,  EMIS 
Datareviews  Series  No.  21,  edited  by  T.  P.  Pearsall  (INSPEC, 
London,  2000). 

61I.  Kudman  and  E.  F.  Steigmeier,  Phys.  Rev.  133,  A1665  (1964). 

62W.  Both,  V.  Gottschalch,  and  G.  Wagner,  Cryst.  Res.  Technol.  21, 
K85  (1986). 


165201-14 


AB  INITIO  THERMAL  TRANSPORT  IN  COMPOUND  . . . 

63P.  V.  Tamarin  and  S.  S.  Shaylt,  Sov.  Phys.  Semicond.  5,  1097 
(1971). 

64G.  Le  Guillou  and  H.  J.  Albany,  Phys.  Rev.  B  5,  2301  (1972). 

65D.  G.  Arasly.  R.  N.  Ragimov,  and  M.  I.  Alief,  Sov.  Phys.  Semicond. 
24,  225  (1990). 

66R.  Bowers,  J.  E.  Bauerle,  and  A.  J.  Cornish,  J.  Appl.  Phys.  30,  1050 
(1959). 

67S.  S.  Shalyt,  P.  V.  Tamarin,  and  V.  S.  Ivleva,  Phys.  Lett.  A  32,  29 
(1970). 

68P.  D.  Maycock,  Solid  State  Commun.  10,  161  (1967). 

69D.  T.  Morelli,  J.  P.  Heremans,  and  G.  A.  Slack,  Phys.  Rev.  B  66, 
195304  (2002). 

70G.  A.  Slack,  R.  A.  Tanzilli,  R.  O.  Pohl,  and  J.  W.  Vandersande, 
J.  Phys.  Chem.  Solids  48,  641  (1987). 

71G.  Nilsson  and  G.  Nelin,  Phys.  Rev.  B  6,  3777  (1972). 

72  A.  Onton,  in  Proceedings  of  the  10th  International  Conference  on 
the  Physics  of  Semiconductors,  edited  by  S.  P.  Keller,  J.  C.  Hensel, 
and  F.  Stern  (US  Atomic  Energy  Commission,  New  York,  1970). 

73T.  Azuhata,  T.  Sota,  and  K.  Suzuki,  J.  Phys.:  Condens.  Matter  7, 
1949  (1995). 


PHYSICAL  REVIEW  B  87,  165201  (2013) 

74D.  Strauch,  B.  Dorner,  and  K.  Karch,  in  Proceedings  of  the 
3rd  International  Conference,  on  Phonon  Physics,  edited  by  S. 
Hunklinger,  W.  Ludwig,  and  G.  Weiss  (World  Scientific,  Singapore, 
1990). 

75P.  H.  Borcherds,  K.  Kune,  G.  F.  Alfreys,  and  R.  L.  Hall,  J.  Phys.  C 
12,  4699  (1979). 

76D.  Strauch  andB.  Dorner,  J.  Phys.:  Condens.  Matter 2,  1457  (1990). 
77M.  K.  Farr,  J.  G.  Traylor,  and  S.  K.  Sinha,  Phys.  Rev.  B  11,  1587 
(1975). 

78P.  H.  Borcherds,  G.  F.  Alfrey,  D.  H.  Saunderson,  and  A.  D.  B. 

Woods,  J.  Phys.  C  8,  2022  (1975). 

79N.  S.  Orlova,  Phys.  Status  Solidi  B  119,  541  (1983). 

80R.  Carles,  N.  Saint-Cricq,  J.  B.  Renucci,  M.  A.  Renucci,  and 
A.  Zwick,  Phys.  Rev.  B  22,  4804  (1980). 

S1D.  L.  Price,  J.  M.  Rowe,  and  R.  M.  Nicklow,  Phys.  Rev.  B  3,  1268 
(1971). 

82D.  W.  Feldman,  J.  H.  Parker,  Jr.,  W.  J.  Choyke,  and  L.  Patrick,  Phys. 
Rev.  173,  787  (1968). 

83M.  Schwoerer-Bohning,  A.  T.  Macrander,  M.  Pabst,  and  P.  Pavone, 
Phys.  Status  Solidi  B  215,  177  (1999). 


165201-15 


